home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_12 / 1012084a < prev    next >
Encoding:
Text File  |  1992-10-12  |  2.2 KB  |  67 lines
  1. /* bc program to calculate Chebyshef economized polynomial
  2. ** for evaluation of sin(x) */
  3. /* use bc -l < listing1.bc to get s() function
  4. ** replace t() with c() to fit cos(x)
  5. ** a(x)/x for atan() */
  6. define t(x) {                           /* sin(x)/x */
  7.         if(x==0)return(1.); /* l'Hopital's rule */
  8.     return (s(x) / x); /* put function to be fit here */ }
  9. define b(x) { /* abs(x) */
  10.     if (x < 0) return (-x);
  11.     return (x); }
  12. define m(x, y) { /* max(x,y) */
  13.     if (x > y) return (x);
  14.     return (y); }
  15. n = 22; /* number of Chebyshef terms - need at least 3 more than
  16. ** are used below */
  17. scale = 40;
  18. /* set precision for dc -must be more than desired result */
  19. p = a(1.) * 4; /* pi */
  20. b = p * .25;
  21. /* upper end of curve fit interval;  b = 7/16 for atan() */
  22. a = -b; /* lower end of interval */
  23. /* chebft adapted from Press Flannery et al */
  24. /* "Numerical Recipes" FORTRAN version */
  25. for (k = 1; k <= n; ++k) {
  26. c[k] = 0;
  27. f[k] = t(c((k - .5) * p / n) * (b - a) * .5 + (b + a) * .5);
  28. }
  29. /* because of symmetry, even c[] are 0 */
  30. for (j = 1; j <= n; j += 2) {
  31.     s = 0;
  32.     q = (j - 1) * p / n;
  33.     for (k = 1; k <= n; ++k) s += c(q * (k - .5)) * f[k];
  34.     /* old style bc requires =+ */
  35.     (c[j] = 2 / n * s); } /* display results */
  36. /* skip even terms, which are 0 */
  37. for (n = 5; n <= 19; n += 2) {
  38.  /* chebpc */
  39.     for (j = 1; j <= n; ++j) d[j] =  e[j] = 0;
  40.     d[1] = c[n];
  41.     for (j = n - 1; j >= 2; --j) {
  42.         for (k = n - j + 1; 1 == 1; --k) {
  43.             s = e[k];
  44.             e[k] = d[k]; 
  45.             if(k==1)break;
  46.             d[k] = d[k - 1] * 2 - s; }
  47.         d[1] = c[j] - s; }
  48.     for (j = n; j >= 2; --j) d[j] = d[j - 1] - e[j];
  49.     d[1] = c[1] * .5 - e[1];
  50.  /* pcshft */
  51.     g = 2 / (b - a);
  52.     for (j = 2; j <= n; ++j) {
  53.         d[j] *= g; /* bc doesn't support multiple assignment */
  54.         g *= 2 / (b - a); }
  55.     for (j = 1; j < n; ++j) {
  56.         h = d[n];
  57.         for (k = n - 1; k >= j; --k) {
  58.             h = d[k] - (a + b) * .5 * h;
  59.             d[k] = h; }
  60.     }
  61.     "Chebyshev Sin fit |x|<Pi/4 coefficients"
  62.         " Maximum Rel Error:"
  63.         m(b(c[n + 2]), b(c[2])) / t(b);
  64.     for (i = 1; i <= n; i += 2)  d[i]; 
  65. }
  66. WRAP_EOF
  67.